What the code does:

Calculates all stats for song sound pressure level (i.e. amplitude) covariates

Next steps

check if magnitude of response to treatments in significant models is explained by body size of lifting power

Load packages

## vector with package names
x <- c( "pbapply", "parallel", "ggplot2", "warbleR", "Rraven", "viridis", "readxl", "rptR", "MCMCglmm", "MuMIn", "corrplot", "lme4", "grid", "gridExtra")

aa <- lapply(x, function(y) {
  
  # check if installed, if not then install 
  if (!y %in% installed.packages()[,"Package"]) 
    install.packages(y) 

  # load package
  try(require(y, character.only = T), silent = T)
})

Functions and parameters

#functions and parameters
knitr::opts_knit$set(root.dir = normalizePath(".."))

knitr::opts_chunk$set(dpi = 50, fig.width = 12, warning = FALSE, message = FALSE) 

# ggplot2 theme
# theme_set(theme_classic(base_size = 20))

cut_path <- "./data/raw/cuts"

treatments <- c("Calibration", "Regular_singing", "Coordination", "After_chase",  
 "Before_playback", "After_playback", "Before_interaction", "After_interaction", "Before_noise", "After_noise")

# iterations for MCMCglmm models
itrns <- 100000

# functions from https://rdrr.io/rforge/rptR/src/R/rpt.mcmcLMM.R
rpt.mcmcLMM <- function(y, groups, CI=0.95, prior=NULL, verbose=FALSE, ...){
    # initial checks
    if(length(y)!= length(groups)) stop("y and group are of unequal length")
    # preparation
    groups <- factor(groups)
    if(is.null(prior)) prior <- list(R=list(V=1,n=10e-2), G=list(G1=list(V=1,n=10e-2)) )
    # point estimation according to model 8 and equation 9
    mod   <- MCMCglmm(y ~ 1, random=~groups, family="gaussian", data=data.frame(y=y,groups=groups), prior=prior, verbose=verbose, ...)
    var.a <- mod$VCV[,"groups"]
    var.e <- mod$VCV[,"units"]
    postR <- var.a / (var.a + var.e)
    # point estimate
    R     <- posterior.mode( postR )
    # credibility interval estimation from paterior distribution
    CI.R    <- coda::HPDinterval(postR,CI)[1,]
    se      <- sd(postR)
    # 'significance test'
    P     <- NA
    res = list(call=match.call(), datatype="Gaussian", method="LMM.MCMC", CI=CI, 
                R=R, CI.R=CI.R, se=se, P=P) 
    # class(res) <- "rpt"
    return(res) 
}


## print Gelman-Rubin convergence statistics, plots traces and autocorrelations
mcmc_diagnostics <- function(rep_mods_list){
  
  for(w in 1:length(rep_mods_list))
{
  mod_name <- names(rep_mods_list)[w]
  
   if(mod_name == "1") mod_name <- "Null" 
  print(paste("model:", mod_name))
  
  Y <- lapply(rep_mods_list[[w]], "[[", "Sol")
  

  ## add global plots and gelman test
  # gelman_diagnostic
  gel_diag <- as.data.frame(gelman.diag(mcmc.list(Y))$psrf)
  
  # add estimate as column
  gel_diag$estimate <- rownames(gel_diag)
  
  # reorder columns
  gel_diag <- gel_diag[, c(3, 1, 2)]

  par(mfrow = c(1, 1))

  # plot table
  grid.newpage()
  grid.draw(tableGrob(gel_diag, rows = NULL, theme=ttheme_default(base_size = 25)))  

  par(mfrow = c(1, 4))

  traceplot(Y, col = adjustcolor(c("yellow","blue", "red"), alpha.f = 0.6))
  
  autocorr.plot(x = Y[[1]], auto.layout = FALSE, lwd =4, col = "red")
} 
  
}

amp <- read.csv("./output/calibrated_amplitude_all_songs.csv")
amp$Treatment[amp$Treatment == "Regular_sining"] <- "Regular_singing"

Repeatability

We estimated repeatability using linear mixed models with MCMC sampling (MCMCglmm() function)

# keep only morning recs, only regular singing and only from 2021
amp_morn <- amp[amp$period == "morning" & amp$year == 2021 & amp$Treatment == "Regular_singing", ]

# keep only males recorded at least twice 
count_sf <- table(amp_morn$ID[!duplicated(amp_morn$org.sound.file)])

amp_morn_rep <- amp_morn[amp_morn$ID %in% names(count_sf)[count_sf > 1], ]

max_quantile <- c(0, 1, seq(10, 90, by = 10)) / 100
outliers <- c(0, 1, 2, 5, 10, 20) / 100

rep_grid <- expand.grid(max_quantile = max_quantile, outliers = outliers, only.low.outliers = c(TRUE, FALSE))

repts <- pblapply(1:nrow(rep_grid), cl = 6, function(x){  
  
  quant_subet <- lapply(unique(amp_morn_rep$org.sound.file), function(y){
    X <- amp_morn_rep[amp_morn_rep$org.sound.file == y, ]
    
    # remove outliers
    outlier_quant <- quantile(X$cal.spl, c(rep_grid$outliers[x], 1- rep_grid$outliers[x]))
  
        if (!rep_grid$only.low.outliers[x])
    X <- X[X$cal.spl >= outlier_quant[1] & X$cal.spl <= outlier_quant[2],] else  X <- X[X$cal.spl >= outlier_quant[1],]
    
    # quantlie for each max quantile
    quant <- quantile(X$cal.spl, probs = 1 - rep_grid$max_quantile[x])
    
    # subset
   X <- X[X$cal.spl >= quant, ]
  })
  
  quant_subet <- do.call(rbind, quant_subet)
  
  # frequentist
    # out <- rpt(cal.spl ~ (1 | ID), grname = "ID", data = quant_subet, datatype = "Gaussian", npermut = 0, nboot = 100)

  # bayesian
    out <- rpt.mcmcLMM(y = quant_subet$cal.spl, groups = quant_subet$ID, nitt = itrns)
    
  out_df <- data.frame(rep_grid$max_quantile[x], rep_grid$outliers[x], rep_grid$only.low.outliers[x], out$R, out$CI.R[1], out$CI.R[2])
  
  return(out_df)
  })

repts_df <- do.call(rbind, repts)

names(repts_df) <- c("max_quantile", "outliers", "only.low.outliers", "repeatability", "lowCI", "hiCI")

write.csv(repts_df, "./output/repeatability_optimization.csv", row.names = FALSE)

Different subsets of data for upper quantiles (x axis)
Pannels show the proportion of outliers removed

repts_df <- read.csv("./output/repeatability_optimization.csv")

pd <- position_dodge(width = 0.1)
ggplot(data = repts_df, aes(x = 1 - max_quantile, y = repeatability, color = only.low.outliers, group = only.low.outliers)) + 
  geom_hline(yintercept = 0.5, col = adjustcolor("red", alpha.f = 0.5)) +
  geom_point(size = 2, position = pd) +
  geom_errorbar(width=.05, aes(ymin = lowCI, ymax = hiCI), position = pd) +
  scale_color_viridis(discrete = TRUE, begin = 0.2, end = 0.8, alpha = 0.7) +
  geom_line(position = pd) + 
  labs(y = "Repeatability",  x = "Upper quantile used") +
  ylim(c(0, 1)) + xlim(c(1, 0)) +
  facet_wrap(~ outliers, scales = "fixed") + 
  theme_classic(base_size = 24)

repts_df$range <- repts_df$hiCI - repts_df$lowCI

# repts_df[order(repts_df$repeatability), c("max_quantile", "outliers", "only.low.outliers", "repeatability", "range")]

repts_df <- repts_df[order(repts_df$repeatability, decreasing = TRUE), ]

kable(head(repts_df))

	max_quantile	outliers	only.low.outliers	repeatability	lowCI	hiCI	range
126	0.3	0.20	FALSE	0.6159727	0.2920327	0.8264643	0.5344317
94	0.4	0.02	FALSE	0.5707252	0.2818027	0.8171755	0.5353727
40	0.5	0.05	TRUE	0.5645533	0.2930169	0.8232646	0.5302477
98	0.8	0.02	FALSE	0.5600182	0.2848066	0.8065382	0.5217316
93	0.3	0.02	FALSE	0.5589800	0.2987828	0.8229877	0.5242049
127	0.4	0.20	FALSE	0.5588000	0.2928387	0.8197835	0.5269448

Excluding lowest values increases repeatability
Excluding outliers doesn’t have a strong effect but removing 0.1 outliers produced the highest repeatability
Removing lower tail outliers or both works similarly (we excluded only low tail outliers to keep more data)

Subsetting the 0.4 upper quartile after excluding 0.2 lower tail outliers:

# compose variable to remove low values and outliers based on repeatabiliy
amp$osf.treat <- paste(amp$org.sound.file, amp$Treatment, sep = "-")

rm_outlier_amp_l <- lapply(unique(amp$osf.treat), function(y){
    X <- amp[amp$osf.treat == y, ]
    
    # remove outliers
    outlier_quant <- quantile(X$cal.spl, c(repts_df$outliers[1], 1))
    X <- X[X$cal.spl >= outlier_quant[1] & X$cal.spl <= outlier_quant[2],]
    
    # quantlie for each max quantile (0.6 was selected due to high repeatability)
    quant <- quantile(X$cal.spl, probs = repts_df$max_quantile[1])
    
    # subset
   X <- X[X$cal.spl >= quant, ]
  })

rm_outlier_amp <- do.call(rbind, rm_outlier_amp_l)

Hypothesis testing analysis

We use bayesian MCMC generalized linear models with individual ID and song type as random factors (random intercepts) except indicated. An intercept-only (null) model is also included in the analyses. A loop is used to run these models. Each model is replicated three times with starting values sampled from a Z-distribution (“start” argument in MCMCglmm()). Diagnostic plots for MCMC model performance are shown at the end of this report.

SPL vs day period (morning / afternoon)

Only for individuals recorded in both periods

Period of the day as predictor: \[SPL \sim + period + (1 | individual) + (1 | song\ type)\]
Null model with no predictor (intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]

# subset data
period_data <- rm_outlier_amp[rm_outlier_amp$Treatment == "Regular_singing" & rm_outlier_amp$year == 2021 & rm_outlier_amp$ID %in% unique(rm_outlier_amp$ID[rm_outlier_amp$period == "afternoon"]), ]

# model SPL by body size
period_formulas <- c("cal.spl ~ 1", "cal.spl ~ period")

period_mods <- pblapply(period_formulas, function(x){
  
  replicate(n = 3, MCMCglmm(as.formula(x), random = ~ ID + songtype, data = period_data, verbose = FALSE, nitt = itrns,  start = list(QUASI = FALSE)), simplify = FALSE)
  
})

names(period_mods) <- gsub("cal.spl ~ ", "", period_formulas)
saveRDS(list(period_formulas = period_formulas, period_mods = period_mods), "./output/mcmcglmm_period_models.RDS")

Model selection

attach(readRDS("./output/mcmcglmm_period_models.RDS"))

mod_list <- lapply(period_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl ~ ", "", period_formulas)

model_selection <- model.sel(mod_list, rank="DIC")

model_selection

	(Intercept)	period	family	df	logLik	DIC	delta	weight
period	98.38209	+	(NA)	5	-6192.009	12389.99	0.00000	1.000000e+00
1	99.28286	NA	(NA)	4	-6223.692	12452.40	62.40826	2.806839e-14

Best model summary

# fixed effects with HPD intervals
best_mod_period <- period_mods[[which(names(period_mods) == row.names(model_selection)[1])]][[1]]

sm <- as.data.frame(summary(best_mod_period)$solutions[, -5])

sm

	post.mean	l-95% CI	u-95% CI	eff.samp
(Intercept)	98.382089	90.2424713	106.17904	4700
periodmorning	1.232936	0.9223919	1.52236	4700

dat_period <- rm_outlier_amp[rm_outlier_amp$Treatment == "Regular_singing" & rm_outlier_amp$ID %in% unique(rm_outlier_amp$ID[rm_outlier_amp$period == "afternoon"]), ]

dat_period$period <- factor(dat_period$period, levels = c("morning", "afternoon"))

ggplot(dat_period, aes(x = as.factor(ID), y = cal.spl, color = period, fill = period)) +
  geom_violin(position =  pd) +
  labs(x = "Individual", y = "Sound pressure level (dB)", fill = "Period") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.6) +
  scale_fill_viridis_d(begin = 0.1, end = 0.8, alpha = 0.6) +
  guides(color = FALSE) +
  theme_classic(base_size = 24)

Effect size posterior distribution

# simplify names
colnames(best_mod_period$Sol) <- c("Intercept", "Morning vs afternoon")

# stack posteriors
Y <- stack(as.data.frame(best_mod_period$Sol[, -1]))
Y$ind <- "Morning vs afternoon"

# plot posteriors
ggplot(Y, aes(y=values, x = ind)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
 geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  labs(y = "Effect size posterior distribution",x = "Predictor") +
  theme_classic(base_size = 24) +
  theme(axis.text.x = element_text(angle = 45,  vjust = 0.9, hjust=1))

Morning songs have significantly higher SPL than afternoon songs within individuals

SPL vs condition

Only for individuals recorded in the morning (afternoon recordings were excluded)
Condition parameters: body size (PC1 from a PCA on total culmen, flattened wing length, body mass and central rectrix) and lifting power

Body size as predictor: \[SPL \sim + body\ size + (1 | individual) + (1 | song\ type)\]
Lifting power as predictor: \[SPL \sim + lifting\ power + (1 | individual) + (1 | song\ type)\]
Lifting power and body size as predictors: \[SPL \sim + lifting\ power + body\ size + (1 | individual) + (1 | song\ type)\]
Interaction of lifting power and body size as predictor: \[SPL \sim + lifting\ power * body\ size + (1 | individual) + (1 | song\ type)\]
Null model with no predictor (intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]

caps <- read_excel("/home/m/Dropbox/LBH data/Additional data files/LBH captures data.xlsx")

# same variables as in paper on spatial memory and body size  
size_vars <- caps[caps$`Bird ID` %in% amp$ID, c("Bird ID", "Total culmen", "Flattened wing length", "Weight", "Central rectriz")]

# replace 419 central rectriz with mean
size_vars$`Central rectriz`[size_vars$`Bird ID` == 419] <- mean(size_vars$`Central rectriz`, na.rm = TRUE)

cor(size_vars[, c("Total culmen", "Flattened wing length", "Weight", "Central rectriz")], use = "pairwise.complete.obs")

##                       Total culmen Flattened wing length     Weight
## Total culmen           1.000000000           -0.13969106 0.14426785
## Flattened wing length -0.139691059            1.00000000 0.08646515
## Weight                 0.144267849            0.08646515 1.00000000
## Central rectriz        0.007536892            0.34107800 0.10235729
##                       Central rectriz
## Total culmen              0.007536892
## Flattened wing length     0.341077996
## Weight                    0.102357292
## Central rectriz           1.000000000

mean_size <- aggregate(. ~ `Bird ID`, data = size_vars, FUN = mean, na.action = "na.omit")

pca <- prcomp(mean_size[, c("Total culmen", "Flattened wing length", "Weight", "Central rectriz")], scale. = TRUE)

summary(pca)

## Importance of components:
##                           PC1    PC2    PC3    PC4
## Standard deviation     1.2488 1.0841 0.8726 0.7098
## Proportion of Variance 0.3899 0.2938 0.1903 0.1259
## Cumulative Proportion  0.3899 0.6837 0.8741 1.0000

mean_size$PC1 <- pca$x[, 1]


rm_outlier_amp$PC1 <- sapply(1:nrow(rm_outlier_amp), function(x) mean_size$PC1[mean_size$`Bird ID` == rm_outlier_amp$ID[x]][1])

rm_outlier_amp$weight <- sapply(1:nrow(rm_outlier_amp), function(x) mean_size$Weight[mean_size$`Bird ID` == rm_outlier_amp$ID[x]][1])

### add weight lifted
lift <- read_excel("/home/m/Dropbox/LBH data/Morphology-condition/other files/Uplift power experiment.xlsx", sheet = "Results")

rm_outlier_amp$lift.weight <- sapply(1:nrow(rm_outlier_amp), function(x){
  
  if(rm_outlier_amp$ID[x] == 0) weight <- NA else {
    
    if (rm_outlier_amp$year[x] == 2019)
      sub_lift <- lift[grep("\\.2019\\.", lift$Video), ] else
              sub_lift <- lift[grep("\\.2021\\.", lift$Video), ]

    
    weights <- sub_lift$Weight[sub_lift$ID == rm_outlier_amp$ID[x]]
   
    if (length(weights) > 0)
    # get the mean of the 2 highest flights
     weight <- mean(sort(weights, decreasing = TRUE)[1:2]) else
    weight <- NA
  }  
  
  return(weight)
})

# subset data
body_size_data <- rm_outlier_amp[rm_outlier_amp$Treatment == "Regular_singing" & rm_outlier_amp$year == 2021 & !is.na(rm_outlier_amp$PC1) & !is.na(rm_outlier_amp$lift.weight) & rm_outlier_amp$period == "morning", ]

# model SPL by body condition
size_formulas <- c("cal.spl ~ 1", "cal.spl ~ PC1", "cal.spl ~ lift.weight", "cal.spl ~ lift.weight + PC1", "cal.spl ~ lift.weight *  PC1")


size_mods <- pblapply(size_formulas, function(x){
  
  replicate(n = 3, MCMCglmm(as.formula(x), random = ~ ID + songtype, data = body_size_data, verbose = FALSE, nitt = itrns,  start = list(QUASI = FALSE)), simplify = FALSE)
  
})

names(size_mods) <- gsub("cal.spl ~ ", "", size_formulas)
saveRDS(list(size_formulas = size_formulas, size_mods = size_mods), "./output/mcmcglmm_size_models.RDS")

Model selection

attach(readRDS("./output/mcmcglmm_size_models.RDS"))

mod_list <- lapply(size_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl ~ ", "", size_formulas)

model_selection <- model.sel(mod_list, rank="DIC")

model_selection

	(Intercept)	PC1	lift.weight	lift.weight:PC1	family	df	logLik	DIC	delta	weight
PC1	101.06722	1.0132676	NA	NA	(NA)	5	-19465.51	38951.92	0.0000000	0.2082655
lift.weight + PC1	85.44995	0.8044757	0.9775365	NA	(NA)	6	-19465.54	38951.92	0.0073858	0.2074978
lift.weight * PC1	77.28863	10.1490394	1.5058787	-0.6219876	(NA)	7	-19465.49	38952.02	0.1043244	0.1976803
lift.weight	84.06902	NA	1.0776599	NA	(NA)	5	-19465.55	38952.06	0.1464912	0.1935562
1	101.35748	NA	NA	NA	(NA)	4	-19465.57	38952.07	0.1522446	0.1930002

Best model summary

None of the models provided a better fit compared to the null model

# fixed effects with HPD intervals
best_mod_size <- size_mods[[which(names(size_mods) ==  "lift.weight + PC1")]][[1]]

sm_size <- as.data.frame(summary(best_mod_size)$solutions[, -5])

sm_size

	post.mean	l-95% CI	u-95% CI	eff.samp
(Intercept)	85.4499479	49.409360	122.753355	4700.00
lift.weight	0.9775365	-1.421323	3.169801	4700.00
PC1	0.8044757	-2.093328	3.909391	4354.76

body_size_data <- rm_outlier_amp[rm_outlier_amp$Treatment == "Regular_singing" & rm_outlier_amp$year == 2021 & !is.na(rm_outlier_amp$PC1) & !is.na(rm_outlier_amp$lift.weight) & rm_outlier_amp$period == "morning", ]

agg_size <- aggregate(cbind(cal.spl, lift.weight, PC1) ~ ID, data = body_size_data, mean)
agg_size$count <- aggregate(PC1 ~ ID, data = body_size_data, length)$PC1

ggplot(agg_size, aes(x = PC1, y = cal.spl)) +
    geom_jitter(position = position_jitter(width = 0.2, height = 0.2), aes(size = count)) +
  labs(x = "Body size (PC1)", y = "Sound pressure level (dB)", color = "Period", size = "Sample size") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24) + 
  # geom_abline(slope = c(ES_continious = sm_size[2, 1], ES_continious_plus_ES_interaction=  sm_size[2, 1] + sm_size[4, 1]), intercept = c(intercept = sm_size[1, 1], intercept_plus_ES_categorical = sm_size[1, 1] + sm_size[3, 1]), color = viridis(2, begin = 0.2, end = 0.8, alpha = 0.6)[2:1], size = 1) +
guides(colour = guide_legend(override.aes = list(size=8)))

ggplot(agg_size, aes(x = PC1, y = cal.spl)) +
    geom_jitter(position = position_jitter(width = 0.2, height = 0.2), aes(size = count)) +
  labs(x = "Body size (PC1)", y = "Sound pressure level (dB)", color = "Period", size = "Sample size") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24) + 
  # geom_abline(slope = c(ES_continious = sm_size[2, 1], ES_continious_plus_ES_interaction=  sm_size[2, 1] + sm_size[4, 1]), intercept = c(intercept = sm_size[1, 1], intercept_plus_ES_categorical = sm_size[1, 1] + sm_size[3, 1]), color = viridis(2, begin = 0.2, end = 0.8, alpha = 0.6)[2:1], size = 1) +
guides(colour = guide_legend(override.aes = list(size=8)))

ggplot(agg_size, aes(x = lift.weight, y = cal.spl)) +
    geom_jitter(position = position_jitter(width = 0.2, height = 0.2), aes(size = count)) +
  labs(x = "Lifiting power (g)", y = "Sound pressure level (dB)", color = "Period", size = "Sample size") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24) + 
  # geom_abline(slope = c(ES_continious = sm_size[2, 1], ES_continious_plus_ES_interaction=  sm_size[2, 1] + sm_size[4, 1]), intercept = c(intercept = sm_size[1, 1], intercept_plus_ES_categorical = sm_size[1, 1] + sm_size[3, 1]), color = viridis(2, begin = 0.2, end = 0.8, alpha = 0.6)[2:1], size = 1) +
guides(colour = guide_legend(override.aes = list(size=8)))

Each dot represents the average for a single individual

Effect size posterior distribution

# simplify names
colnames(best_mod_size$Sol) <- c("Intercept", "lifted weight","Body size (PC1)")

# stack posteriors
Y <- stack(as.data.frame(best_mod_size$Sol[, -1]))

# plot posteriors
ggplot(Y, aes(y=values, x = ind)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
 geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  labs(y = "Effect size posterior distribution", x = "Predictor") +
  theme_classic(base_size = 24) +
  theme(axis.text.x = element_text(angle = 45,  vjust = 1, hjust=1))

Song PSL is not affected by condition (either body size or lifting power)

SPL vs playback

Using data from 2019 and 2021
Recordings from 2019 were not calibrated and were center around SPL mean from 2021
Comparing before playback vs after playback (during regular singing)

Playback as predictor and random slope: \[SPL \sim + body\ size + (playback | individual) + (1 | song\ type)\]
Playback as predictor and random slope: \[SPL \sim + body\ size + (playback | individual) + (1 | song\ type)\]
Playback as predictor and only random intercept: \[SPL \sim + playback + (1 | individual) + (1 | song\ type)\]
Null model with no predictor (random intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]

playback_dat <- rm_outlier_amp[rm_outlier_amp$Treatment %in% c("Before_playback", "After_playback"), ]

# label 2019 recordings
playback_dat$ID <- ifelse(grepl("2019", playback_dat$org.sound.file), paste(playback_dat$ID, "(2019)"), playback_dat$ID)

# remove indiviudal recorded in 2019 and 2021
playback_dat <- playback_dat[playback_dat$ID != "397 (2019)", ]

random_slope_mod <- lmer(cal.spl ~ Treatment  + (Treatment|ID), data = playback_dat, REML = FALSE)

interaction_size_mod <- lmer(cal.spl ~ PC1 + Treatment + PC1:Treatment + (1|ID), data = playback_dat, REML = FALSE)

interaction_lift_mod <- lmer(cal.spl ~ PC1 + lift.weight + PC1:lift.weight + (1|ID), data = playback_dat, REML = FALSE)

random_intercept_only_mod <- lmer(cal.spl ~ Treatment  + (1|ID) + (1|songtype), data = playback_dat, REML = FALSE)

null_mod <- lmer(cal.spl ~ 1  + (1|ID) + (1|songtype), data = playback_dat, REML = FALSE)

d1 <- (-2*logLik(random_intercept_only_mod)) + (2*logLik(random_slope_mod))

pval_rand_slope_vs_rand_intercept <- 0.5*pchisq(d1[1], 0, lower.tail=FALSE) + 0.5*pchisq(d1[1], 1, lower.tail=FALSE)

d2 <- (-2*logLik(null_mod)) + (2*logLik(random_intercept_only_mod))

pval_rand_intercept_vs_null <- 0.5*pchisq(d2[1], 0, lower.tail=FALSE) + 0.5*pchisq(d2[1], 1, lower.tail=FALSE)

d3 <- (-2*logLik(null_mod)) + (2*logLik(random_slope_mod))

pval_rand_slope_vs_null <- 0.5*pchisq(d3[1], 0, lower.tail=FALSE) + 0.5*pchisq(d3[1], 1, lower.tail=FALSE)

d4 <- (-2*logLik(interaction_size_mod)) + (2*logLik(random_slope_mod))

pval_rand_slope_vs_interaction_size <- 0.5*pchisq(d4[1], 0, lower.tail=FALSE) + 0.5*pchisq(d4[1], 1, lower.tail=FALSE)

d5 <- (-2*logLik(interaction_lift_mod)) + (2*logLik(random_slope_mod))

pval_rand_slope_vs_interaction_lift <- 0.5*pchisq(d5[1], 0, lower.tail=FALSE) + 0.5*pchisq(d5[1], 1, lower.tail=FALSE)

agg_playback <- aggregate(cal.spl ~ Treatment + ID, data = playback_dat, mean)

agg_playback$sd <- aggregate(cal.spl ~ Treatment + ID, data = playback_dat, sd)$cal.spl

agg_playback$Treatment <- factor(gsub("_playback", "", agg_playback$Treatment), levels = c("Before", "After"))

ggplot(agg_playback, aes(x = Treatment, y = cal.spl, color = Treatment)) +
  geom_point(size = 2, show.legend = FALSE) +
   geom_errorbar(width=.05, aes(ymin = cal.spl - sd, ymax = cal.spl + sd), show.legend = FALSE) +
  labs(x = "Playback treatment", y = "Sound pressure level (dB)", color = "Period") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 1) +
  facet_wrap(~ ID, nrow = 3, scales = "free_y") +
  theme_classic(base_size = 24)

The random slope model is significantly better than the random intercept model (p = < 0.0001), the interaction with body size model (p = < 0.0001), the interaction with lifting power model (p = < 0.0001) and the null model (p = < 0.0001)

The playback stimuli elicits a significant change in SPL but the direction of the response varies between individuals
Direction of the response doesn’t seem to be affected by body size

SPL vs coordination

Random slope models, using data from 2019 and 2021
Only for individuals in which coordination was observed when recorded
SPL of songs sang during coordinated singing were compared to regular singing songs

Coordination as predictor and random slope: \[SPL \sim + coordination + (coordination | individual) + (1 | song\ type)\] 1.Coordination as predictor and random intercept only: \[SPL \sim + coordination + (1 | individual) + (1 | song\ type)\]
Null model with no predictor (random intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]

# fixed effects with HPD intervals
coord_data <- rm_outlier_amp[rm_outlier_amp$org.sound.file %in% rm_outlier_amp$org.sound.file[rm_outlier_amp$Treatment == "Coordination"], ]

coord_data <- coord_data[coord_data$Treatment %in%  c("Regular_singing","Coordination"), ]

random_slope_mod <- lmer(cal.spl ~ Treatment  + (Treatment|ID) + (1|songtype), data = coord_data, REML = FALSE)

random_intercept_only_mod <- lmer(cal.spl ~ Treatment + (1|ID) + (1|songtype), data = coord_data, REML = FALSE)

null_mod <- lmer(cal.spl ~ 1  + (1|ID) + (1|songtype), data = coord_data, REML = FALSE)

d1 <- (-2*logLik(random_intercept_only_mod)) + (2*logLik(random_slope_mod))

# 0.5*pchisq(d1[1], 0, lower.tail=FALSE) + 0.5*pchisq(d1[1], 1, lower.tail=FALSE)


d2 <- (-2*logLik(null_mod)) + (2*logLik(random_intercept_only_mod))

# 0.5*pchisq(d2[1], 0, lower.tail=FALSE) + 0.5*pchisq(d2[1], 1, lower.tail=FALSE)

d3 <- (-2*logLik(null_mod)) + (2*logLik(random_slope_mod))

# 0.5*pchisq(d3[1], 0, lower.tail=FALSE) + 0.5*pchisq(d3[1], 1, lower.tail=FALSE)

# sm_playback <- as.data.frame(fixef(random_slope_mod), confint(random_slope_mod))
# 
# sm_playback

summary(random_slope_mod)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: cal.spl ~ Treatment + (Treatment | ID) + (1 | songtype)
##    Data: coord_data
## 
##      AIC      BIC   logLik deviance df.resid 
##  21588.1  21631.8 -10787.1  21574.1     3768 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9918 -0.1442  0.0093  0.1648  3.8007 
## 
## Random effects:
##  Groups   Name                     Variance  Std.Dev. Corr 
##  ID       (Intercept)              7.513e+01 8.667822      
##           TreatmentRegular_singing 3.656e+00 1.912156 -0.26
##  songtype (Intercept)              4.264e-06 0.002065      
##  Residual                          1.730e+01 4.158752      
## Number of obs: 3775, groups:  ID, 11; songtype, 9
## 
## Fixed effects:
##                          Estimate Std. Error t value
## (Intercept)              101.7875     2.6173  38.890
## TreatmentRegular_singing  -0.2087     0.6293  -0.332
## 
## Correlation of Fixed Effects:
##             (Intr)
## TrtmntRglr_ -0.251

agg_coord <- aggregate(cal.spl ~ Treatment + ID, data = coord_data, mean)

agg_coord$sd <- aggregate(cal.spl ~ Treatment + ID, data = coord_data, sd)$cal.spl

agg_coord$Treatment <- factor(ifelse(agg_coord$Treatment == "Coordination", "Coordinated", "Solo"), levels = c("Solo", "Coordinated"))

ggplot(agg_coord, aes(x = Treatment, y = cal.spl, color = Treatment)) +
  geom_point(size = 2, show.legend = FALSE) +
   geom_errorbar(width=.05, aes(ymin = cal.spl - sd, ymax = cal.spl + sd), show.legend = FALSE) +
 
  labs(x = "coord treatment", y = "Sound pressure level (dB)", color = "Period") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 1) +
  facet_wrap(~ ID, nrow = 2, scales = "free_y") +
  theme_classic(base_size = 24)

The random slope model was not significantly better than the random intercept model (p = < 0.0001) nor the null model (p = < 0.0001)

Vocal coordination is not associated with detectable changes in SPL

SPL vs aggresive interaction

2 types of interactions: perch interaction and chases
SPL of songs before (i.e. regular singing) and after interactions were compared
Interactions encoded in 2 ways:
- 1. factor clumping perch interactions and chases in a single category (2 levels: before and after interaction)
- 1. Interaction type (chase of perch interaction)

Interaction as predictor: \[SPL \sim + interaction + (1 | individual) + (1 | song\ type)\]
Interaction and interaction type as predictors: \[SPL \sim + interaction + interaction\ type + (1 | individual) + (1 | song\ type)\]
Statistical interaction between “Interaction” and “interaction type” as predictor: \[SPL \sim + interaction * interaction\ type + (1 | individual) + (1 | song\ type)\]
Null model with no predictor (intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]

# subset data
dat_chase <- rm_outlier_amp[rm_outlier_amp$org.sound.file %in% rm_outlier_amp$org.sound.file[rm_outlier_amp$Treatment == "After_chase"], ]

dat_chase <- dat_chase[dat_chase$Treatment %in% c("Regular_singing", "After_chase"),]


dat_chase$inter_treatment <- ifelse(grepl("After", dat_chase$Treatment), "After", "Before")

dat_chase$inter_type <- "Chase"

dat_interaction <- rm_outlier_amp[rm_outlier_amp$Treatment %in% c("After_interaction", "Before_interaction"), ]


dat_interaction$inter_treatment <- ifelse(grepl("After", dat_interaction$Treatment), "After", "Before")
dat_interaction$inter_type <- "Perch_interaction"

# fix range of non-calibrated recs
dat_interaction$cal.spl[dat_interaction$year < 2021] <- dat_interaction$cal.spl[dat_interaction$year < 2021] + mean(dat_interaction$cal.spl[dat_interaction$year == 2021]) - mean(dat_interaction$cal.spl[dat_interaction$year < 2021])

dat_aggresive_inter <- rbind(dat_interaction, dat_chase)

# fix low value one
dat_aggresive_inter$cal.spl[dat_aggresive_inter$sound.files == dat_aggresive_inter$sound.files[which.min(dat_aggresive_inter$cal.spl)]] <- dat_aggresive_inter$cal.spl[dat_aggresive_inter$sound.files == dat_aggresive_inter$sound.files[which.min(dat_aggresive_inter$cal.spl)]]  +  mean(dat_aggresive_inter$cal.spl[dat_aggresive_inter$sound.files != dat_aggresive_inter$sound.files[which.min(dat_aggresive_inter$cal.spl)]]) - mean(dat_aggresive_inter$cal.spl[dat_aggresive_inter$sound.files == dat_aggresive_inter$sound.files[which.min(dat_aggresive_inter$cal.spl)]])

# table(dat_aggresive_inter$ID, dat_aggresive_inter$inter_treatment, dat_aggresive_inter$inter_type)

# model SPL by interaction
interaction_formulas <- c("cal.spl ~ 1", "cal.spl ~ inter_treatment", "cal.spl ~ inter_treatment + inter_type", "cal.spl ~ inter_treatment * inter_type")

inter_mods <- pblapply(interaction_formulas, function(x){
  
  replicate(n = 3, MCMCglmm(as.formula(x), random = ~ ID + songtype, data = dat_aggresive_inter, verbose = FALSE, nitt = itrns,  start = list(QUASI = FALSE)), simplify = FALSE)
  
})

names(inter_mods) <- gsub("cal.spl ~ ", "", interaction_formulas)

saveRDS(list(interaction_formulas = interaction_formulas, inter_mods = inter_mods), "./output/mcmcglmm_interaction_models.RDS")

Model selection

attach(readRDS("./output/mcmcglmm_interaction_models.RDS"))

mod_list <- lapply(inter_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl ~ ", "", interaction_formulas)

model_selection <- model.sel(mod_list, rank="DIC")

model_selection

	(Intercept)	inter_treatment	inter_type	inter_treatment:inter_type	family	df	logLik	DIC	delta	weight
inter_treatment * inter_type	105.8312	+	+	+	(NA)	7	-2926.655	5874.963	0.00000	1.000000e+00
inter_treatment + inter_type	104.1285	+	+	NA	(NA)	6	-2960.277	5941.127	66.16386	4.292404e-15
inter_treatment	102.7296	+	NA	NA	(NA)	5	-2973.693	5966.872	91.90835	1.102439e-20
1	101.4128	NA	NA	NA	(NA)	4	-3051.179	6120.967	246.00364	3.810561e-54

Best model summary

# fixed effects with HPD intervals
best_mod_interaction <- inter_mods[[which(names(inter_mods) == row.names(model_selection)[1])]][[1]]

sm_interaction <- as.data.frame(summary(best_mod_interaction)$solutions[, -5])

sm_interaction

	post.mean	l-95% CI	u-95% CI	eff.samp
(Intercept)	105.831196	103.633244	108.344906	1993.982
inter_treatmentBefore	-4.252822	-4.765444	-3.725284	4700.000
inter_typePerch_interaction	-3.970494	-4.775029	-3.131514	4700.000
inter_treatmentBefore:inter_typePerch_interaction	3.042986	2.351168	3.764623	4700.000

# subset data
dat_aggresive_inter$inter_treatment <- factor(dat_aggresive_inter$inter_treatment, levels = c("Before", "After"))

dat_aggresive_inter$year <- factor(dat_aggresive_inter$year, levels = c("2013", "2014", "2019", "2021"))

dat_aggresive_inter$inter_type <- gsub("_", " ", dat_aggresive_inter$inter_type)
  
dat_aggresive_inter$inter_type <- factor(dat_aggresive_inter$inter_type, levels = c("Perch interaction", "Chase"))

agg_aggresive_inter <- aggregate(cal.spl ~ ID + inter_type + inter_treatment + year, data = dat_aggresive_inter, mean)

agg_aggresive_inter$ID <- as.character(agg_aggresive_inter$ID)

ggplot(dat_aggresive_inter, aes(x = inter_treatment, y = cal.spl)) +
  geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  geom_point(size = 4, data = agg_aggresive_inter, aes(group = ID, color = year)) +
  geom_line(data = agg_aggresive_inter, aes(group = ID, color = year)) +
  labs(x = "Context", y = "Sound pressure level (dB)") +
  facet_wrap(~ inter_type) +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24)

Effect size posterior distribution

# simplify names
colnames(best_mod_interaction$Sol) <- c("Intercept", "After vs before interaction", "Chase vs perch_interaction", "After vs before chase")

# stack posteriors
Y <- stack(as.data.frame(best_mod_interaction$Sol[, -1]))

# plot posteriors
ggplot(Y, aes(y=values, x = ind)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
 geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  labs(y = "Effect size posterior distribution",x = "Predictor") +
  theme_classic(base_size = 24) +
  theme(axis.text.x = element_text(angle = 45,  vjust = 0.9, hjust=1))

Male-male aggressive interactions result in a significant SPL increase
Chases have a stronger effect, increasing SPL in ~ 4 dB compare to regular singing before the chase

SPL and song structure

Song structure measured as spectrographic consistency and gap duration

Song spectrographic consistency

spectrographic consistency: mean spectrographic cross-correlation of all songs in a singing bout
Spectrographic consistency was measured for all songs within contexts (before & after interaction, clumping chase and perch interactions)
Spectrographic consistency was compared between the contexts
Recordings from interactions were used as they showed the highest increase in SPL

interaction as predictor and spectrographic consistency as response: \[consistency \sim + interaction + (1 | individual) + (1 | song\ type)\]
interaction and SPL as predictors and spectrographic consistency as response: \[consistency \sim + interaction + SPL + (1 | individual) + (1 | song\ type)\]
interaction between “interaction” and SPL as predictor and spectrographic consistency as response: \[consistency \sim + interaction * SPL + (1 | individual) + (1 | song\ type)\]
Null model with no predictor (intercept only): \[consistency \sim + 1 + (1 | individual) + (1 | song\ type)\]

agg_amp_inter <- aggregate(cal.spl ~ ID + inter_type, data = dat_aggresive_inter[dat_aggresive_inter$inter_treatment == "After",], mean)  

agg_amp_inter$amp.change <- agg_amp_inter$cal.spl -  
aggregate(cal.spl ~ ID + inter_type, data = dat_aggresive_inter[dat_aggresive_inter$inter_treatment == "Before",], mean)$cal.spl  

# sort by change
agg_amp_inter <- agg_amp_inter[order(agg_amp_inter$amp.change, decreasing = TRUE), ]

# subset data
dat_song_consistency <- dat_aggresive_inter[dat_aggresive_inter$ID %in% agg_amp_inter$ID[agg_amp_inter$amp.change >= 0.5], ]

dat_song_consistency$ID.treatment <- paste(dat_song_consistency$ID, dat_song_consistency$inter_treatment, sep ="-")

dat_song_consistency_l <- lapply(unique(dat_song_consistency$ID.treatment), function(x){
  
  X <- dat_song_consistency[dat_song_consistency$ID.treatment == x, ]
  
  print(x)
  if(nrow(X) > 1){
        # measure cross-correlation
    xcrr <- (cross_correlation(X, path = "./data/raw/cuts", pb = FALSE, wl = 150))  
    
    # measure gaps
    gps <- gaps(X, pb = FALSE)$gaps
    gps <- gps[gps < 1]
    mean.gaps <- mean(gps, na.rm = TRUE)
    
    # mean xcorr
    if (!is.na(xcrr[1])){
    xcrr <- mean(xcrr[lower.tri(xcrr)])
  }
  
  out_df <- data.frame(sound.files = x, ID = X$ID[1], Treatment = X$Treatment[1], xcrr = xcrr,  mean.gaps = mean.gaps, mean.spl = mean(X$cal.spl), songtype = X$songtype[1], treatment = X$inter_treatment[1])
  } else 
  out_df <- data.frame(sound.files = x, ID = X$ID[1], Treatment = X$Treatment[1], xcrr = NA,  mean.gaps = NA, mean.spl = mean(X$cal.spl), songtype = X$songtype[1], treatment = X$inter_treatment[1])
  
  return(out_df)
  
})

dat_song_consistency <- do.call(rbind, dat_song_consistency_l)

dat_song_consistency <- dat_song_consistency[!is.na(dat_song_consistency$xcrr), ]

saveRDS(dat_song_consistency, "./output/dat_song_spectrographic_consistency_and_gaps.RDS")

dat_song_consistency <- readRDS("./output/dat_song_spectrographic_consistency_and_gaps.RDS")


# model SPL by body size
song_consistency_formulas <- c("xcrr ~ 1", "xcrr ~ treatment", "xcrr ~ treatment + mean.spl", "xcrr ~ treatment * mean.spl" , "xcrr ~mean.spl")

song_consistency_mods <- pblapply(song_consistency_formulas, function(x){
  
  replicate(n = 3, MCMCglmm(as.formula(x), random = ~ ID + songtype, data = dat_song_consistency, verbose = FALSE, nitt = itrns,  start = list(QUASI = FALSE)), simplify = FALSE)
  
})

names(song_consistency_mods) <- gsub("cal.spl ~ ", "", song_consistency_formulas)

saveRDS(list(song_consistency_formulas = song_consistency_formulas, song_consistency_mods = song_consistency_mods), "./output/mcmcglmm_song_consistency_models.RDS")

Model selection

attach(readRDS("./output/mcmcglmm_song_consistency_models.RDS"))

mod_list <- lapply(song_consistency_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl ~ ", "", song_consistency_formulas)

model_selection <- model.sel(mod_list, rank="DIC")

model_selection

	(Intercept)	treatment	mean.spl	mean.spl:treatment	family	df	logLik	DIC	delta	weight
xcrr ~ treatment	0.7383678	+	NA	NA	(NA)	5	34.73809	-62.44242	0.000000	0.62783916
xcrr ~ treatment + mean.spl	0.2875213	+	0.0043802	NA	(NA)	6	32.91965	-59.39476	3.047666	0.13679057
xcrr ~ 1	0.7546914	NA	NA	NA	(NA)	4	32.56895	-59.33922	3.103208	0.13304401
xcrr ~mean.spl	0.5292037	NA	0.0022100	NA	(NA)	5	31.75161	-57.88083	4.561592	0.06416694
xcrr ~ treatment * mean.spl	0.2022626	+	0.0051972	+	(NA)	7	31.84366	-56.84140	5.601029	0.03815931

Best model summary

# summary
best_mod_song_consistency <- song_consistency_mods[[which(names(song_consistency_mods) == row.names(model_selection)[1])]][[1]]

sm_song_consistency <- as.data.frame(summary(best_mod_song_consistency)$solutions[, -5])

sm_song_consistency

	post.mean	l-95% CI	u-95% CI	eff.samp
(Intercept)	0.7383678	0.6867783	0.7959804	4689.774
treatmentBefore	0.0323433	-0.0122364	0.0807822	9700.000

dat_song_consistency <- readRDS("./output/dat_song_spectrographic_consistency_and_gaps.RDS")

# subset data
dat_song_consistency$ID <- as.character(dat_song_consistency$ID)

ggplot(dat_song_consistency, aes(x = mean.spl, y = xcrr)) +
  geom_point(size = 4, color = viridis(5, alpha = 0.5)[3]) +
  labs(y = "Mean cross-correlation", x = "Sound pressure level (dB)") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24)

Effect size posterior distribution

# stack posteriors
Y <- stack(as.data.frame(best_mod_song_consistency$Sol[, -1]))
Y$ind <- "mean.spl"

# plot posteriors
ggplot(Y, aes(y=values, x = ind)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
 geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  labs(y = "Effect size posterior distribution",x = "Predictor") +
  theme_classic(base_size = 24) +
  theme(axis.text.x = element_text(angle = 45,  vjust = 0.9, hjust=1))

Aggressive interaction but not SPL affects the spectrographic consistency of songs

Song gap duration

Gaps: silent time intervals between songs

interaction as predictor and gap duration as response: \[gap\ duration \sim + interaction + (1 | individual) + (1 | song\ type)\]
interaction and SPL as predictors and gap duration as response: \[gap\ duration \sim + interaction + SPL + (1 | individual) + (1 | song\ type)\]
interaction between “interaction” and SPL as predictor and gap duration as response: \[gap\ duration \sim + interaction * SPL + (1 | individual) + (1 | song\ type)\]
Null model with no predictor (intercept only): \[gap\ duration \sim + 1 + (1 | individual) + (1 | song\ type)\]

dat_gaps <- readRDS("./output/dat_song_spectrographic_consistency_and_gaps.RDS")

# model SPL by body size
gap_formulas <- c("mean.gaps ~ 1", "mean.gaps ~ inter_treatment", "mean.gaps ~ inter_treatment + mean.spl", "mean.gaps ~ inter_treatment * mean.spl" , "mean.gaps ~ mean.spl")

gap_mods <- pblapply(gap_formulas, function(x){
  
  replicate(n = 3, MCMCglmm(as.formula(x), random = ~ ID + songtype, data = dat_gaps, verbose = FALSE, nitt = itrns,  start = list(QUASI = FALSE)), simplify = FALSE)
  
})

names(gap_mods) <- gsub("cal.spl ~ ", "", gap_formulas)

saveRDS(list(gap_formulas = gap_formulas, gap_mods = gap_mods), "./output/mcmcglmm_gap_models.RDS")

Model selection

attach(readRDS("./output/mcmcglmm_gap_models.RDS"))

mod_list <- lapply(gap_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl ~ ", "", gap_formulas)

model_selection <- model.sel(mod_list, rank="DIC")

model_selection

	(Intercept)	inter_treatment	mean.spl	inter_treatment:mean.spl	family	df	logLik	DIC	delta	weight
mean.gaps ~ inter_treatment * mean.spl	-0.0181886	+	0.0067237	+	(NA)	7	69.40951	-127.7444	0.0000000	0.3578651599
mean.gaps ~ inter_treatment	0.6695200	+	NA	NA	(NA)	5	68.40630	-127.6362	0.1081806	0.3390223244
mean.gaps ~ inter_treatment + mean.spl	0.2087143	+	0.0044737	NA	(NA)	6	68.92134	-127.3879	0.3564624	0.2994432988
mean.gaps ~ 1	0.7031442	NA	NA	NA	(NA)	4	63.05805	-118.0621	9.6822948	0.0028264143
mean.gaps ~ mean.spl	0.6270569	NA	0.0007432	NA	(NA)	5	62.37306	-115.6420	12.1023576	0.0008428026

Best model summary

# summary
best_mod_gap <- gap_mods[[which(names(gap_mods) == row.names(model_selection)[1])]][[1]]

sm_gap <- as.data.frame(summary(best_mod_gap)$solutions[, -5])

sm_gap

	post.mean	l-95% CI	u-95% CI	eff.samp
(Intercept)	-0.0181886	-0.7200648	0.7339871	2975.767
inter_treatmentBefore	0.5859397	-0.1818819	1.3476811	9700.000
mean.spl	0.0067237	-0.0000158	0.0140674	3457.288
inter_treatmentBefore:mean.spl	-0.0051698	-0.0129531	0.0022992	9700.000

dat_gaps <- readRDS("./output/dat_song_spectrographic_consistency_and_gaps.RDS")

# subset data
dat_gaps$ID <- as.character(dat_gaps$ID)

ggplot(dat_gaps, aes(x = mean.gaps, y = xcrr)) +
  geom_point(size = 4, color = viridis(5, alpha = 0.5)[3]) +
  labs(y = "Mean gap duration (s)", x = "Sound pressure level (dB)") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24)

Effect size posterior distribution

# stack posteriors
Y <- stack(as.data.frame(best_mod_gap$Sol[, -1]))

# plot posteriors
ggplot(Y, aes(y=values, x = ind)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
 geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  labs(y = "Effect size posterior distribution",x = "Predictor") +
  theme_classic(base_size = 24) +
  theme(axis.text.x = element_text(angle = 45,  vjust = 0.9, hjust=1))

# zoom in
Y <- stack(as.data.frame(best_mod_gap$Sol[, -c(1, 2)]))
ggplot(Y, aes(y=values, x = ind)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) +
 geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25, begin = 0.4)) +
  labs(y = "Effect size posterior distribution",x = "Predictor") +
  theme_classic(base_size = 24) +
  theme(axis.text.x = element_text(angle = 45,  vjust = 0.9, hjust=1))

No detectable effect on gap duration

Diagnostic plots for MCMCglmm models

Include Gelman/Rubin’s convergence diagnostic, MCMC chain trace (all tree replicates in a single plot: yellow, blue and red colors) and autocorrelation plots:

SPL and period of the day

mcmc_diagnostics(period_mods)

## [1] "model: Null"

## [1] "model: period"

SPL and condition

mcmc_diagnostics(size_mods)

## [1] "model: Null"

## [1] "model: PC1"

## [1] "model: lift.weight"

## [1] "model: lift.weight + PC1"

## [1] "model: lift.weight *  PC1"

SPL and interactions

mcmc_diagnostics(inter_mods)

## [1] "model: Null"

## [1] "model: inter_treatment"

## [1] "model: inter_treatment + inter_type"

## [1] "model: inter_treatment * inter_type"

SPL and spectrographic consistency

mcmc_diagnostics(song_consistency_mods)

## [1] "model: xcrr ~ 1"

## [1] "model: xcrr ~ treatment"

## [1] "model: xcrr ~ treatment + mean.spl"

## [1] "model: xcrr ~ treatment * mean.spl"

## [1] "model: xcrr ~mean.spl"

SPL and gap duration

mcmc_diagnostics(gap_mods)

## [1] "model: mean.gaps ~ 1"

## [1] "model: mean.gaps ~ inter_treatment"

## [1] "model: mean.gaps ~ inter_treatment + mean.spl"

## [1] "model: mean.gaps ~ inter_treatment * mean.spl"

## [1] "model: mean.gaps ~ mean.spl"

Session information

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=pt_PT.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_PT.UTF-8    
##  [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_PT.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] gridExtra_2.3      lme4_1.1-26        corrplot_0.84      MuMIn_1.43.17     
##  [5] MCMCglmm_2.29      ape_5.4-1          coda_0.19-4        Matrix_1.2-18     
##  [9] rptR_0.9.22        readxl_1.3.1       viridis_0.5.1      viridisLite_0.3.0 
## [13] Rraven_1.0.10      warbleR_1.1.27     NatureSounds_1.0.3 knitr_1.31        
## [17] seewave_2.1.6      tuneR_1.3.3        ggplot2_3.3.2      pbapply_1.4-3     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6        lattice_0.20-41   corpcor_1.6.9     fftw_1.0-6       
##  [5] digest_0.6.27     R6_2.5.0          cellranger_1.1.0  signal_0.7-6     
##  [9] stats4_4.0.3      evaluate_0.14     highr_0.8         pillar_1.4.6     
## [13] rlang_0.4.10      cubature_2.0.4.1  minqa_1.2.4       nloptr_1.2.2.2   
## [17] rmarkdown_2.4     labeling_0.3      splines_4.0.3     statmod_1.4.35   
## [21] stringr_1.4.0     RCurl_1.98-1.3    munsell_0.5.0     proxy_0.4-25     
## [25] compiler_4.0.3    xfun_0.22         pkgconfig_2.0.3   htmltools_0.5.1.1
## [29] tidyselect_1.1.0  tibble_3.0.4      tensorA_0.36.1    dtw_1.22-3       
## [33] crayon_1.4.1      dplyr_1.0.2       withr_2.3.0       MASS_7.3-53      
## [37] bitops_1.0-6      nlme_3.1-149      gtable_0.3.0      lifecycle_0.2.0  
## [41] magrittr_2.0.1    scales_1.1.1      stringi_1.5.3     farver_2.0.3     
## [45] ellipsis_0.3.1    generics_0.1.0    vctrs_0.3.6       boot_1.3-25      
## [49] rjson_0.2.20      tools_4.0.3       glue_1.4.2        purrr_0.3.4      
## [53] yaml_2.2.1        colorspace_1.4-1

Sound pressure level stats

Long-billed hermit song amplitude

Marcelo Araya-Salas, PhD & Melanie Talavera

19-04-2021

What the code does:

Next steps

Load packages

Functions and parameters

Repeatability

Hypothesis testing analysis

SPL vs day period (morning / afternoon)

Model selection

Best model summary

SPL vs condition

Model selection

Best model summary

SPL vs playback

SPL vs coordination

SPL vs aggresive interaction

Model selection

Best model summary

SPL and song structure

Song spectrographic consistency

Model selection

Best model summary

Song gap duration

Model selection

Best model summary

Diagnostic plots for MCMCglmm models

SPL and period of the day

SPL and condition

SPL and interactions

SPL and spectrographic consistency

SPL and gap duration